% Code for ``Towards a framework for a new research ecosystem"
% by Andrea Modena and Roberto Savona. This version: 30/05/2025
% This code simulates the performance of the R&D investment portfolio in
% Figure 3.

% Technical parameters
T =10;              % Total time
N = T*250;           % Number of time steps
dt = T/N;           % Time step size
t = linspace(0, T, N+1);  % Time vector
M = 5000;             % Number of Brownian motion paths


% Model parameters
sig=0.15; % risky asset volatility
mu=0.085; % risk asset expected return
r=0.046; % risk-free rate
gamma=2; % risk aversion
rho=0.02; % return target
lam=0.04; % innovation intensity
G=0.1; % innovation return
th0=0.015; % time 0 R&D spending
PM = poissrnd(lam*dt, M, N); % innovation events

% Rootfinder

% Define the function
f = @(y) (exp(y*T) - 0.3 - y/th0); % finds the growth rate of R&D cost
% Find the root (exponent of the growth rate of phi)
delta =0.12;

% Initial guess (you can tweak this depending on expected solution)
x0 = 0.2;



% Simulate Brownian increments
dW = sqrt(dt) * randn(M, N);  % M paths × N increments
W = [zeros(M,1), cumsum(dW, 2)];  % Start at zero, then accumulate

E=zeros(M,N+1);  
h=zeros(M,N+1);  

E(:,1)=1; % Sets fund initial value


% Simulates fund value paths
for i=1:M
    for j=2:N
        h(i,j)=exp(rho*(T))*exp(-r*(T-t(j)))-lam*G*(1-exp(-r*(T-t(j))))-exp(delta*t(j))*(th0/(delta-r))*(1-exp(-(r-delta)*(T-t(j))));
        TH_S(i,j)=(E(i,j-1)-h(i,j))*(mu-r)/(gamma*sig^2);
        E(i,j)=E(i,j-1)+E(i,j-1)*r*dt+TH_S(i,j)*dt*(mu-r)-th0*exp(delta*t(j))*dt...
            +TH_S(i,j)*sig*dW(i,j)+G*PM(i,j);
    end
end

E=E(:,1:end-1);
t=t(1:end-1);

figure;
% Plot paths
plot(t, E', 'Color', [0.2 0.6 1 0.1], 'LineWidth', 0.1);
hold on

% Plots mean and percentiles
plot(t,mean(E)','b--','LineWidth', 1.5);
plot(t,prctile(E, 95),'b:','LineWidth', 1.5);
plot(t,prctile(E, 5),'b:','LineWidth', 1.5);
grid on
xlabel('Time (years)', 'Interpreter', 'latex')
ylabel('Fund value', 'Interpreter', 'latex')
